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Error Analysis of p- version discontinuous Galerkin method 
for heat transfer in built-up structures 

H. Kaneko and K. S. Bey 
Abstract 

The purpose of this paper is to provide an error analysis for the p- version of the discon- 
tinuous Galerkin finite element method for heat transfer in built-up structures. As a special 
case of the results in this paper, a theoretical error estimate for the numerical experiments 
recently conducted by James Tomey [7] is obtained. 
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1 Introduction. 

The purposes of this paper are to report the state of the art information on time discretization 
techniques in the discontinuous Galerkin method for parabolic problems (this section) and to 
establish error analysis for p- version of the finite element method for such problems (section 
2). Also a discussion of various time discretization techniques are included (section 3). The 
discontinuous Galerkin method is applied to the following standard model problem of parabolic 
type: 

Find u such that 

u t (x, t) — A u(x,t) = f{x,t), x € ft,t £ R+, 

u(x,t) = 0, x € dft,t £ R+, (1-1) 

u(x,0) — uo(z), x £ ft, 

where 0 is a dosed and bounded set in R ! 3 with boundary dft, R+ = (0, oo), Au = d 2 u/dx 2 -f 
d 2 /dy 2 + d 2 u/dz 2 , u t = du/dt , and the functions / and uo are given data. In this paper, region 
ft is assumed to be a thin body in R 3 , such as a panel on the wing or fuselage of an aerospace 
vehicle, p-version of the finite element method is considered in all directions including time 
variable. Because of the special characteristics of the region, it is assumed that, through the 
thickness, only one element is taken. This allows us to avoid construction of finite elements 
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that are too thin to violate the quasi uniform condition (see (1.11) below). The recent paper [5] 
addresses the similar issue under the framework of the 'modified’ hp- finite element scheme. For 
a Banach space A and / = (0, T) indicating a time interval, we denote by 1^(1; X), 1 < p < oo, 
and H k (I;X ), 0 < k e R, the Lebesgue and Sobolev spaces. Also P P (I;X) denotes the set 
of all polynomials of degree < p with coefficients in A, i.e., q(t) e P P {I\X) if and only if 
q(t) = Yfj=o x jV for some xj € X and tel . Let Tj be a partition of I into M (I) subintervals 
{In — (tn-i,t n )}n=P ' We set k n = t n - t n - 1 . Denote by and u~ the right and the left limits 
of u at t n respectively. We set [u]„ = u+ -u~,n= 1, . . . , M(I) - 1. For each time interval J n , 
an approximation order p„ > 0 is assigned and they are stored in the vector p = The 

semidiscrete space is then given by 


V p (T r , X) = {u: I —+ X: u|/ n € P p ”(I n ; A), 1 < n < M(J)}. 


If p is a constant vector, i.e., p n = p for all 1 < n < M(I ), then we write V p (Tj; A). The 
number of degrees of freedom of the time discretization will be denoted by NRDOFfV^fTV; A)) 
and it is, of course, equal to J2n=i(Pn + !)• The semidiscrete solution U e P Pn (I n ; A) of the 
problem (1.1), if U is already determined on /&, 1 < k < n — 1, is found by solving the following 
problem: 

Find U e P Pn (/ n ; A) such that 

[ m,V) 2 + ( V U, V V) 2 }dt + (U+_ 1 ,V+_ 1 ) 2 = [ (g,V)x'x X dt + (U~_ 1 ,V r ^_ 1 ) 2 (1.2) 

•'in J I n 

for all V e P Pn (I n ; A) and = uq- 

Here we assume that L 2 (Q) is densely embedded in a Banach space A. The following theorem 
and its subsequent corollary are reported in [10]. Theorem describes the error estimate for the 
semidiscrete solution explicitly in terms of the time steps, the approximation orders and the 
local regularities of the solution. 


Theorem 1.1 Let u be the solution of (1.1) and U the semidiscrete solution in V p {Tj\ A). 
Assume that u\j n e fP n+1 (J n ;A) for 1 < n < M(I) and s n , 1 < n < M(I), nonnegative 
integers. Then 






Hlr.-X) 


for any 0 < s* < min(p n ,s„). 
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In the p- version of the discontinuous Galerkin finite element method, a temporal partition Tj is 
fixed and convergence is obtained by p n — > oc. For p n = P for each n = 1, . . . , M(I), and for a 
smooth solution u, we obtain the following: 

Corollary 1.2 Letp n = p /or each n = 1, . . . M(I) and k max = ma x{k n }. Let u E H So+l (I; X) 7 
for a nonnegative integer so, be the exact solution of (LI) and U E V p (Tj\X) the semidiscrete 
solution . Then 

u 2{min(5 0 ,p)+l} 

U U*) < 

where C depends on so, but is independent o//r max and p. 

Remark 1.1: As pointed out in [10], Corollary 1.2 shows that for smooth solutions where 
so is large, it is better to increase p than to reduce k max at a fixed, often low p. Since N = 
NRDOF(V^(I/; X)) ~ p, we see that for p-version of the finite element method, 

llu - Ui LHI , X) < CN - ». (1.3) 

Using the standard approximation theory for analytic functions (1.3) reduces to 

II u-U\\ L 2 {I . x) <Ce-<*, (1-4) 

for some b > 0 independent of p. If the solution is not smooth in time, it is still possible to 
approximate it in exponential orders by a /ip-finite element method which combines a certain 
geometric partition with the semidiscrete space V p (Tj;X) where p is linearly distributed, (see 
[10] for detail). Using the h- finite element method with non-uniform graded time partitions, 
such non-smooth solutions can be approximated in an algebraically optimal order, but not 
exponentially, using different approaches, (see, e.g., [5], [10]). The standard p-finite element 
method does not perform well in this context. Therefore, since we aim to establish the p- version 
of the finite element method for (1.1), we assume for the remaining of this paper that solutions 
are smooth in time. This assumption allows us to establish the order of approximation in time 
variable that is compatible with the approximation orders in the spatial variables. 

Now we consider the problem of discretizing the space Q. For simplicity, we assume that 



where u? is divided into M(h) number of triangular elements, each denoted by K l . Let K ^ 
denote the master triangular element defined by 

= {(tv) € fl 2 : 0 < i) < (1 + 0V5 -1 < i < 0 or 

0<7?< (l-C)VS 0<£< 1}. 

Let S P (K ^) denote the space of polynomials of degree < p on -i.e., 

S P (K^) = span{i l rf: i,j = 0, 1, . . . ,p; i + j < p}. 

First, the shape functions for the master element are formed. To accomplish this, the 
barycentric coordinates are introduced via 

Ai = (1 - £ - *?/V3)/2, A 2 = (l+e-r ? /v / 3)/2, A 3 = rf/y/Z. 

Aj’s form a partition of unity and A* is identically equal to one at a vertex of v and vanishes 
on the opposite side of The hierarchical shape functions on consists of internal as well 
as external functions. The normalized antiderivatives of the Legendre polynomials are defined 

by 

= j^Pi(t)dt, i = 1,2,3,... 

Now, the external shape functions consist of 3 nodal shape functions 

Ni(tv) = i= 1,2,3, 

and 3(p — 1) side shape functions 1 — 1, . . . ,p — 1, j = 1,2,3. The index j indicates 

one of three sides of K^. Noting that ^(±1) = 0, 

Mv) = j(l ~V 2 )Vi(v), i= 1,2,3,... (1.4) 

where <^(77) is a polynomial of degree i - 1. For instance, ^1(77) = — y/6, <p 2 ( 77 ) = — y/lQr] and 

^3(77) = —2(1 — 5 p 2 ), etc. The side shape functions are constructed as follows: 

A^(f , V) — (A 3 - A 2 ) 

Af 1 ^) = A 3 Ai<pi(Ai - A 3 ), i — 1, . . . ,p — 1, (l- 5 ) 

N\ Z] (tv) -A 1 A 2 ^(A 2 -A 1 ). 

From (1.4) and (1.5), there are 3 -f 3(p - 1) shape functions. As dim(S p (K^ v )) — 

the remaining ) basis elements are constructed in terms of internal shape functions. 
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Clearly, nontrivial internal shape functions on K^ v exists only if p > 3. For p — 3, the bubble 
function on defined below serves as the internal function; 

V) = Ai A2A3 = ^((1 - ^) 2 - <S 2 ). ( 1 . 6 ) 

Moreover, the collection I P (K '^) of higher-order internal shape functions can be constructed 
from 

I p (Kt,v) = {**<., v: v € S p - 3 (AT C ,„)} - {&*<,„} ® S P ' 3 (K^), p > 3. 

Let Th, h > 0, be a triangulation of u. let x ~ Q l x {L \ , L 2 , L 3 ) and y = Q^(L 1} L 2 , L 3 ) be the 

element mappings of the standard triangle K^ v to the /th triangular element K l E Th, e.g. } the 

linear mapping onto K l with vertices {(x[, yj)}? = i, 

3 3 

Q 1 x{Li,L 2 ,L 3 ) = Q l y (L u L 2 ,L 3 ) = 

i = 1 i=l 

The space of all polynomials of degree < p on K l is denoted by S p (K l ) and its basis can be 
formed from the shape functions of S P (K^^ V ) described above by transforming them under Q l x 
and Q l y . The finite element space S p ' l (u,Th) is now defined. For uj, p > 0 and k > 0, 

S p ' k {u,T h ) = {ue H k (u>):u\ K € S P (K), K € T h }. (1.7) 

Assume that a triangulation {Th}. h > 0, of a; consists of {K l h and that hi = diam(K l h ), 
for / — 1, . . . ,M(h). 

In the z-variable for through the thickness approximation, the local variable r is defined in 
the reference element [—1,1] and T is mapped onto the reference element by Q z , i.e., 

r*Q,([-i,i]), z = q 2 (t). 

Clearly, Q z is a linear function defined by 

2 = Q z (t) = 1(1 - r)(-~) + 1(1 + r)^, r e [-1, 1] 

The Jacobian of Q z is constant 

dz _ d 
dr 2 

In this paper, the basis functions of P p ([— 1, 1]) are taken to be the one-dimensional hierarchical 
shape functions. See [8] for a complete discussion of the basis elements used in the p and hp- finite 
element methods. 
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For example, in approximating an element in H l [— 1,1], with l = 0, ipi(r) = P*_i(r), 1 < 
i < p - hi, where P 7 _i is the Legendre polynomial of degree i — 1, form the hierarchical basis 
functions. With l — 1, the external (xp\ and tfo) and internal (xpi, i > 3) shape functions are 
defined by 

Mr) = Mr) = ^ (lg) 

M T ) = (^) 1/2 /-i Pi-2(t)dt, 3<i<p+l. 

Note that xpiS form an orthogonal family with respect to the energy inner product (•, -)e, 

{A, Me = J 1 M t )M t ) dt = j x Pi(t)Pj(t)dt = Sij. 

Also note that the internal shape functions satisfy 

^{(±1) =0, for 3 < i < p -f 1. 

For the case l — 2 and p > 3, the four nodal shape functions and the remaining p - 3 internal 
shape functions given by 

Mr) = 3(1 - r) 2 ( 1 + r), M T ) = 3(1 - r ) 2 { 2 + r) 

A{r) = —3(1 + r) 2 (l - r), Mr) = 3(1 + r) 2 (2 - t) (1.9) 

A(r) = i^) 1 ' 2 X!j(r - rj)Pi- 3 {v)dv, i = 5, . . . ,p + 1. 

In this case, the internal shape functions satisfy 

dP xb ■ 

— -(±1) =0, for 5 < i < p + 1 and j = 0, 1. 

ar-? 

For example, using the shape functions in (1.8), any element u € H l [— 1, 1] can be approximated 
by u p e P p ([— 1, 1]), in the form 

M t ) = “(-!) + + 5Z a iA( T )- (110) 

z z i= 3 

For approximating the solutions of parabolic problems with the homogeneous Dirichlet boundary 
condition, the first two terms will be dropped, as u(— 1) = u(l) = 0. A sequence of triangulations 
{Th}h>o is called the quasiuniform mesh if 

fora "' ,>0 - < in > 
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with h = maxKer h diam(K), and for some 7 > 0. P p ( T) denotes the space of all polynomials of 
degree < p defined on T. 

The following is proved by Babuska, Szabo and Katz in [1]. See also [2] by Babuska and Suri 
on a related discussion. Here fto denotes a bounded polygonal domain in R 2 . 

Theorem 1.3 Let u € H k (Sla). Then there exists a sequence z p e P p (Uq), p — 1,2,..., such 
that, for any 0 < l < k, 

\\u - Zp||i,n 0 ^ Cp~( k ~ 1 ^ ||u||fc f noi 
where C is independent of u and p. 

The parameters k and l are not necessarily integral. Their proof relies heavily on the approxi- 
mation power of the trigonometric polynomials. 

With l = 0 in Theorem 1.3 and using the usual duality argument, the results in Theorem 
1.3 are further improved by Babuska and Suri in [2] (theorem 2.9), (see also a series of papers 
by Gui and Babuska [9]), to the hp- finite element setting as follows: 

Theorem 1.4 Let T k be a quasiuniform partition o/fto- Then for k > 1, u G H k (Qo), 

v h 2 (no) ^ Ch u p~ k \\u 
veSP*(u),T h ) v ' 

where v = min (k,p T 1)- 

Note that Corollary 1.2 can be derived also from Theorem 1.4, and they establish the same 
result in terms of time variable approximation and spatial variable approximation, respectively. 
The corresponding error estimate in the || • ||//*(n 0 ) a * so available in [2]. Now assume that Ho 

in Theorems 1.3 and 1.4 is ft = u; x [— and consider the problem of approximating elements 
in H k (u; x [—3, |]) by the tensor product space 5 p 1 ,/c (cj, Th) 0 P^— f]. Using Theorems 1.3 

and 1.4, the following is proved in [5]: 

Theorem 1.5 Let u 6 H k {u x |]). Then there exists u* £ S Pl,fc (u/,Tk) 0 P P2 ([-^, 5]), 

11 “ - «1 = °( hl,1 PT k + fc ), 

where Ui = niiiiffc.p, + 1) for i = 1,2 and h = max/cg t h diam(K), with Th a triangulation of w. 
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A straightforward extension of Theorem 1.5 provides the following which describes a total error 
estimate for approximating functions in H k ( u> x [~|,|] x /) by elements from S Pl ' k (oJ,Th) ® 

P K [-j4]®P p HI). 

Theorem 1.6 Let u € H k (u x [-|, 5] x /). Then there exists u € S pi ’ k (u>, Th) ® J 3 ** 2 ([ — 5, 5]) ® 

P P3 {I), 

H“ - «ii i2( a,x[-f ,f]x/) = ' fc ), 

ic/iere i/* — min {h,pi + 1), i = 1,2,3 and h = max^eT,, cfiam(K) 7 uhih Th a triangulation of lo 
and kjjiax = m.axj<j2<jV &n- 

Proof: Take X in V P (I; X) used in Corollary 1.2 to be H k (ujx[-^, ^]). Let u* € V P3 (7; H k (u>x 
^])) be the element approximating u. From Corollary 1 . 2 , 

11 “ ~ = 0(kmaxP3 *)• 

Approximate u * ! ; d d- t by u by Theorem 1.5. Then, 

1 UJX [— 2 J J 

11“ ~ “II L 2 (wx[-|,|]x/) = 11“ ~ “*lli 2 (uix[-f ,f]x/) + ll“* “ “Hz, 2 (wx[-§,f]xJ) 

= c 1 /i^ P r fc + c 2 ^ P 2 fc + czk% ax pi\ 

where C\, C 2 and C3 are constants independent of h , fc max , Pi, P 2 and P3. 


Using (1.4), we obtain 


Corollary 1. 7 Let u € tf°°(u;x[— |]xJ). T/ien there exists u* € S Pl,fc (u>, 2^)0 P P2 ([— f, 5])® 

M u u II L 2 (^x[-|,f]x/) - oe * 

where C and b are independent ofp\,p 2 and p^. 


2 Version of Discontinuous Galerkin Finite Element Method 

In this section, p- version of discontinuous Galerkin finite element method for the parabolic 
problem (1.1) is described. The main goal here is to provide an error analysis for the p-finite 
element method using the results from Section 1. The semidiscrete approximation equation 
(1.2) is now upgraded to a fully discretized equation below. It is assumed as in Theorem 1.5 
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and Corollary 1-6 that the degree vectors in space, through-the-thickness and time are assumed 
to be pi , p 2 and P 3 respectively. Define 

l^(Pi,P2,P3) = { V:R+ S*’ k ®pn[-^]:V |/ B € P w (/„), n = 1, 2, . . . ,N}, (2.1) 

where 

P3 J i 

P P3 (In) = {V(t) = J2 V&: Vi € S”" k <g> P^l- 2 ]}- 

i— 0 

Then the fully discretized discontinuous Galerkin method can be described as follows: 

Find U e W r ( pi,P2,P3 ) such that for n = 1, 2, • * • , 

[ {(Ut,V) + (vU,vV)}dt + ([U} n „ 1 ,V+_ 1 )= [ ( f,v)dt fcraDVgH^**), (2.2) 

*'/n •'/n 


We consider in (1.1) only the case of isotropic materials along with Dirichlet boundary 
conditions. However, extensions to anisotropic materials as well as mixed boundary conditions 
where Neumann boundary conditions are incorporated are possible and the present analysis 
carries over to these cases. The thesis by Tomey [7] treats transversely anisotropic materials as 
well as isotropic materials along with a mixed boundary condition. 

The notation of [7] is followed. The solution u is approximated over K l x (— |,|)x/ n using 
the outer tensor products: 

«l - (4>®ip®0) T a n = X T a n , (2.3) 

which is 

pi p2 P3 

K l x(-%4)xl n (x,y)4>p(z)d' r {t)alp r 

a=0 / 3=0 7=0 

Equation (2.2) becomes 


E «eT h Ik‘ Jr Ik x( t ) T + VX(Vx) T dt] + x(Ci)x(Ci) T dzdu, ■ a" 

= Efc'er h Ik‘ Jr fl n fxdtdzdw + f K , f r x(*"_ 1 )x(t~_ 1 ) T - dzdu ■ a"' 1 . 

Equation (2.4) can be written in the following matrix form [7] 
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(2.5) 


E (C JC «+K jr , + M jr ,)-a"= E (H^+M^.a"- 1 ), 

K l eT h K l eT h 

where C K i, K K i and are, respectively, the element capacitance, conductance (stiffness) 
and mass matrices, whereas H K i is the element load vector. Thesis by Tomey [7] describes in 
detail as to how these matrices should be assembled. In the next section, we consider different 
bases in time variable, delineating an advantage of each choice of basis elements. The following 
theorem can be proved by minor modifications to the proof of theorem 1.1, [3] and by making 
use of Theorem 1.6. 

Theorem 2.1 Let u € H k (u> x [— ~) x I). Suppose that there is a constant 7 such that the 

time steps k n satisfy k n < 7fc n +i, n = 1,...,AT — 1 and let U n denote the solution of (2.2) 
approximating u att n . Then there is a constant C depending only on 7 and a constant (3, where 
PK > (3hx and px is the diameter of the circle inscribed m K for all K € 7k, such that for 
n — 1, 2, ... ,1V, 

ii« - un K (u * { -i ,f]x/) ^ c ( l + lQ g ^) 1/2 {cw‘ P r fc + * + ^3 k% axP 3 k }, 

where Vi = min {k,pi + 1), i = 1,2,3; C\, C 2 and C3 are constants independent of h, d 7 

Pi, P 2 and p^. 

Similarly, Corollary 1.7 implies the following 

Corollary 2.2 Let u 6 x [— |] x I). Suppose that there is a constant 7 such that the 

time steps k n satisfy k n < 7fc ra +i, n — 1 , . . . , N — 1 and let U 71 denote the solution of (2.2) 
approximating u att n . Then there is a constant C depending only on 7 and a constant (3, where 
PK > 0hx and px is the diameter of the circle inscribed in K for all K € T^, such that for 
n = 1 , 2 ,..., AT, 

where C and b are independent of p\,p 2 and p^. 

Remark 2.1: Maintaining throughout computation a certain accuracy of numerical solution 
obtained from (2.1) is always desirable and Theorem 2.1 gives an insight to the following adaptive 
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scheme. Suppose that a tolerance of S is required. Then, as h. d and k max are known a priori , 
the minimum approximation degrees 7 requirement in spatial, through the thickness and time 
variable are obtained from 


h^ k <6, d^ k <6, k% axP ; k <6. 


(2.6) 


Moreover, for the case of infinitely differentiable u, with p = min(pi . P2-P:i) where p,'s satisfy 
inequalities in (2.6), Corollary 2.2 implies that 


\\u-U n \ 


|]x/) - 


< Ce~ bp < 6. 


Remark 2.2: Thus far, the case for the constant degree vectors was considered. Nonconstant 
degree vectors can also be incorporated easily from Theorem 1.1. For instance, nonconstant 
degree vector p = (p n ) can be derived from the inequality in Theorem 1.1. Note that the bound 
given in Theorem 1.1 combines all errors in time approximation over M(I) intervals. Thus, for 
each n = 1, 2 , . . . , M(J), the error in disrectization in time variable over I n is given by 


11“ - t'felfeX, s C (!)*<'"*» (2.6) 

for any 0 < s* < min(p„, s„). Since ~ P n 2Sn as p„ —* oo, (2.6) becomes 


u-U\\ Li{In , x) <Ck‘ n " +l p-^ +l l 


As k n ' s are known, construct the degree vector p by requiring each component p n to satisfy 


^n n+ Vn (5n+1) < S - 


A construction of nonconstant degree vector corresponding to a triangulation T ^ for the region 
lj is similar. 


3 Discretizations in Time Variable. 

In this section, effects of the use of different basis elements to approximate the solution in time 
variable are considered. The solution u is approximated over K l x (~|, |) x I n using the outer 
tensor products: 


To better illustrate the choice of basis elements for time variable, we write a semidiscrete solution 
U n eP Pn (In,X) by 

Pn 

U n — ^ u j,n@j<n* (3.2) 

j=0 

Here Uj, n £ X are unknown coefficients to be determined. V n is defined similarly. For con- 
venience, we let I = I n and drop the time step index n . Substituting U n of (3.2) and the 
corresponding V n into (1.2), we obtain the following: 

Find {uj}^ =0 c X such that 

52 {[/ O' j 6idt + 6f(t n ~i)0?(t n -. 1 ))(u j ,v i )x+ j 9j6idt(Vuj,S7vi)x} 
i,j= 0 Jl Jl 


= f (s, ViOi)x-*xdt + {U n _ u ViOf (tn-i)} for all {ui}?=o c x - 

i = 0 JI 


(3.3) 


Transforming the integrals into the reference interval [—1,1] under F 1 and letting (see [11]) 


Aij = efadt + 0+(-l)<9+(-l), Bij = j l 9j9id&, 

fl(vi) = {J {go F)6idt, Vi)x, fi(vi) = (t/“_ 1 fi i (-l),u i )x, 

equation (3.3) becomes: 

Find the coefficient C X such that for all {vj}^ C X 


52 Aijiiij^Jx 4- ^BijiViijiVvJx = 52 ^fi(vi) + fi(vi). (3.4) 

0 Z i= 0 ^ 


Here k = k n . The strong form of equation (3.4) is 
p £ 

]T + - BijAuj = -/} + /?, i = 0, 1, . . - p, (3.5) 

i=o 

where o F)§idt and if = U“_ 1 0i(— 1). Hence, in order to execute the p-version of 

the finite element method, the following integrals must be computed for assembly of the global 
matrix. 




Odt . 


The set of the canonical polynomials #*,+i(£) = t v were used in [7] and the components of the 
matrices which represent Aij and Bij were computed exactly. We consider two other alternatives 
for 6. 
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First, we consider 0„(t) = : Let U\ In = Since U n ~ l + = ${?(*), 

the advantage of this basis is that the term (£7 n_1 » + , t ,n “ 1,+ ) in (2.1) is simplified. Also, 


and 


f^ n n/dQ \T 


0 

0 

0 


v 

p+i 

p+2 


o -L- -2_ 

p+l p+2 


-E. . 

2p J 


rtn 

/ 1 


0$ T dt = k n 


1 I 

2 3 

I I 

3 4 


1 

P+1 

1 

p+2 


_L J_ J_ 

L p+l p+2 p+3 


2p+l J 


1 



Odt ~ k n 


1 

2 

l 

3 


A 


As stated in [11], the ideal choice for 9{ would be the one under which the matrices A and B 
diagonalize simultaneously. If the diagonalizations of A and B are possible, then equation (3.4) 
decouples into p - f 1 independent equations, reducing the size of computation. The canonical 
basis and the basis just considered generate the full matrices for A and for B. In order to select 
basis functions in time variable which takes into account of the structures of A and B, we now 
consider the Legendre polynomials for 6 . 

Second: The translated Legendre polynomials for 6 v (t) are considered. This is essentially 
the same as the normalized Legendre polynomials used in [11]. We extend the discussion in 
[11] by exhibiting the general form for A and state its characteristics. The orthogonal nature 
of the Legendre polynomials guarantees the matrix B to be diagonal. Hence it remains to 
analyze the matrix A. In [11], it is stated that \.,this (diagonalization) seems not to be possible 
with time shape functions in R, but numerical experiments show that A ... is diagonalizable 
in C at least for 0 < p < 100’. We are unable to provide a mathematical proof of their 
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statement at this point, but our investigation thus far indicates that A is non-defective (see, 
[4]), -i.e., the algebraic multiplicity and the geometric multiplicity of each eigenvalue are the 
same. Hence, A is diagonalizable over the complex field C. In this paper, we propose to use the 
real Schur decomposition of the matrix A, rather than the diagonalization technique in [11], to 
establish a solution scheme for (3.4). The real Schur decomposition leads to a modified backward 
substitution scheme. The method is mathematically justifiable for any degree p. Moreover, the 
current method avoids complex number arithmetic which was necessary in [11], and thus the 
cost of computation is approximately the same as [11]. 

The advantage of this choice as basis elements in time variable lies in the formations of 
Jtn-i It^_i &dt as well as / t * i 0(^) T dt, all of which are of banded structures as seen 

below. Recall that Pi(x) denote the Legendre polynomial of degree i defined over [—1,1]. Let 

x n (t) = t — j — t - tn ^ tn - 1 and L?(t) = P;(x n (£)) for each n = 1,2, . . . ,M(7) and i = 0,1, 

Obviously, 

- f " U>(t\L n At\dt = Si* for each i, j = 0, 1 , . . . , . 

t n - t„- 1 7t„_! * ' ' J ' ' 


Thus, 


and 


Also 



ee T dt = 



\ J 


ftn 

/ Odt = 
Jtn - 1 


0 



0 k n 0 k n 0 k n * * ■ 

0 0 k n 0 k n 0 

0 0 0 k n 0 k n - • 

0 0 0 0 k n 0 • • ■ 


The formation of the last matrix Q(^) T dt requires some tedious but straightforward calcu- 
lations which are not so obvious. Thus, we include them below. Clearly, it is sufficient to derive 
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the following: 


f 1 , 

/_/<*> dt ‘ 


0 2 0 2 0 2 
0 0 2 0 2 0 
0 0 0 2 0 2 
0 0 0 0 2 0 


(3.6) 


For the Legendre polynomials of the first kind, Pi(x), i > 0, x € [—1, 1], first note that 
P 2i (- 1) = -P 2 i(l) = 1, P 2 i+i (— 1) = -1, and P 2i+ i(l) = 1 for i > 0. 

To derive (3.6), we must show that 

J Pi(t)Pj(t)dt = 0, for all i and j with i > j > 0, 


and 


f'mru.*w~ 0 . J = 2, for all i > 0 and l > 0. 


(3.7) 


(3.8) 


Equations (3.7) and (3.8) are verified by induction. As Po(t) = 1, (3-7) is verified with i = 0. 
Also 

j 1 P 0 (t)^,(t)dt = Po(t)P 2 i(t)l-i - f Pa(t)P 2 i(t)dt = 0, for all l > 0. 
and similarly 

f Po(t)P 2 i+i(t)dt = P 0 (t)P 2 i + i(t)\U - £ P^(t)P 2l+1 (t)dt = 2, for all l > 0. 

Now assume that (3.7) and (3.8) are satisfied for some i = i* — 1 and for all j such that 
0 and for all l > 0. First, for the diagonal element, 

f Pi*(t)P[.(t)dt = Pi*(t)P i *(f)|i 1 - £ PAt)P!.(t)dt 

implies that P t *(t)P'*(t)dt = 0. Also, for i* > j > 0, 


f Pi’(t)Pj(t)dt = Pi*(t)Pj(t) ILj. - £ P'.(t)P j (t)dt = 


2-2 = 0 for i* — j odd, 

0-0 = 0 for i* — j even 
where the inductive assumption was used in computing the second integral. This shows (3.7) for 
all i and i > j > 0. It remains to prove that P (t)P/» +2i (t)dt = 0 and P \+ (t)P/* +2 j +1 (t)cft = 
2 for all l > 0. The case for l = 0 is done. For / > 0, 


J Pi* (t)Pi*+2l(t)dt ~ Pi*{t)Pi*+2l{t)\-l ~~ J ^ Pi*(t)Pi*+2l{t)dt. 
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The last integral is 0 by (3.6) and the term P i *(t)P i * + 2 i{t)\l 1 is 0 regardless of i* even or odd. 
Hence Pi*(t)P(*+ 2l (t)dt — 0. f* x Pi*(t)P/. +2/+1 (^)dt — 2 is similar. Thus, (3.6) is verified. 
Now, from (3.6), it is clear that 

1 1 1 111 

-1 1 1 111 

1-1 1 111 

-1 1-1 111 

1-1 1-111 

Formula (3.9) provides a general construction method for the assembly of the p-finite element 
matrix with any order p when the Legendre polynomials are used in time variable. For p = 5, 
with normalizing factors incorporated, (3.9) can be used to derive the matrix A in [11], (eq. 
(4.11)). Two observations on the matrix A p+1 are as follows: 

Theorem 3.1 Let A p+ i denote the (p + 1) x (p + 1) leading principle matrix of (3.6) for each 
p = 0, 1, — Then A p+ i is invertible for each p. 

Proof: This follows immediately by noting that det{A p +\) = 2 P . 

The matrix A p+ i is not positive definite for some p. However: 

Theorem 3.2 The matrix A p+ i is non-negative definite for all n — 0,1, 

As in (3.2), denote the semidiscrete solution U as well as V in (1.2) as 

u = Y J u j e h v = '£v j e j (3.io) 

j=0 j — 0 

where the subscript n is dropped. Let $i = Li, the transformed Legendre polynomial of degree i 
and $i — Pi , the Legendre polynomial of degree i defined over [-1, 1] in equations (3.3) and (3.4). 
Also, denote a basis for the finite element space S pljk (co, Th) ® |] for X by { Sj with 

D = dim(S pl ' k (u,T h ) 0 PP 2 [— | , |]). The trial and test functions uj and in the semi-discrete 
system (3.3) above are further approximated by 

D D 

u i ~ J2 Uj Sl (x), v!fs k (x). (3.11) 

l=i k= l 
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Substituting (3.11) into (3.4), the fully discrete p-finite element system can be written as 
(ref. [11]), for the unknown coefficient vectors Uj = (uj, w?, . . . , u®) T G R D , 


iooM + |S • 

AopM 


Uq 

_ k 

ft 

4- 

'ft ' 

ApoM 

■ ■ AppM + 1 5 


Up 

~~ 2 

■ tc? 

.ft . 


where 


M = {(sj,a*)x}ft=i, 5 = {(Vs,, Vs k ) x }? k=1 , 

and 

f} = (fj(s 1 )JJ(s 2 ),...Jj(s D )) T 
fj = (/ j 2 ( si ),/ j 2 ( s 2 ) 1 -- - ,fj{sD)) T - 

Note that the use of the translated Legendre polynomial served well because By — S(j in 
(3.5). Equation (3.5) can be written in vector form as 

k 

A u + -[6 ij }Au = F, (3.13) 

where F = (| fi^g o x n )6 0 di + t^floGl), • • • , f fli(9 ° x n )d p di + 1)) T '. In [11], A 

is diagonalized and equation (3.13) is solved for w = Q T u from 

Q t AQw + ^ A w = Q t F. (3.14) 

It is reported in [11] that, as the result of their numerical experiments, the matrix A p+ \ is 
diagonalizable up to its order p = 100. Subsequently, equation (3.13) is decoupled via (3.14) 
into p + 1 independent scalar equations, each of which requires complex arithmetic to solve. A 
new approach which uses the real Schur decomposition is now presented. The new approach 
does not require the complex arithmetic. 


Theorem 3.3 (Real Schur Decomposition , [4] p-219) If A € R nxn , then there exists an orthog- 
onal matrix Q 6 R nxn suc h that 


q t aq = 


#11 #12 • • • #lm 
0 #22 ' ' • #2m 


0 0 


(3.15) 
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where each Ra is either a lxl matrix or a 2x2 matrix having complex conjugate eigenvalues 
of A. 


By Theorem 3.1, it is guaranteed that no 1 x 1 matrix in Schur decomposition for A is 0. 
Also, it is worth noting that the Schur decomposition is an orthogonal similarity transformation 
and thus avoid the computation of Q -1 as required in the diagonalization process done in 
[11]. Equation (3.14) is then solved by the backward substitution using block matrices. More 
specifically, w p -\ and w p are found from 


(RmmMA-S) 


Wp- 1 


w n 


F v - 


p - i 


or from 


{RmmM T — S)uJp — F p , 


(3.16) 


(3.17) 


according to Rmm being 2 x 2 or 1 x 1 matrix respectively. Here, we used the standard convention 


that if Rmn 


a b 
c d 


, then RmmM = 


. Computation proceeds to find w p -2 as 


aM bM 
cM dM 

well as w p - 3 if Rm- i m -x is 2 x 2. Namely, assuming that Rmm is 2 x 2, if Rm- im -i is a scalar, 
then w p _ 2 is found by solving 





k 

Rm— 1 m— T —S 

Wp— 2 — Rp—2 [Rm— 1 mAT] 

W p — 1 
Wp 


If Rm-im-i is 2 x 2, then w p s and w p 2 are found from 


k 

Rjn 1 m — 1 AI T —S 


Wp- 3 

= 

Fp- 3 

- [Rm-lmM] 

Wp— 1 

l 


4 s * 

1 

to 

i 


F P _ 2 


Wp 


The case for 1 x 1 Rm m is similar. It is important to recall that the computation thus described 
can be completed because of Theorem 3.1. 


4 Start-up Singularities: 

In this final section, we make some comments on the start-up singularities normally associated 
with the parabolic problems. The regularity assumptions in Theorem 2.1 and Corollary 2.2 were 
taken so that the current fully p-finite element method could provide numerical solutions where 
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the discretization error associated in time can be made consistent with the discretization error 
associated in space. However, as indicated earlier, time singularities arise due to various types 
of incompatible initial data. To capture such singularities, hp-version of finite element method 
must be considered. In [10], a nonuniform time discretization is determined by considering the 
conditions on / as well as the initial function uo in (1.1). More specifically, the function / is 
assumed to be piecewise analytic as a function on [0, T] with values in H , i.e., 

||/ (i) ||tf < Cd l T(l + 1), t € [0, T], l e No 

with constants C and d independent of l and t. Also, uo is assumed to be in He = (H,X)e t 2 , 
0 < 0 < 1, where (H,X)e ,2 is a space between H and X determined by the if -met hod of inter- 
polation, (see [8]). An h - version of the discontinuous Galerkin finite element method developed 
in [5] establishes a nonuniform time discretization scheme which is based upon the behavior of 
l|u (1) Wllx- The direct inspection into the smoothness of ||u^^(t)||x in determining time parti- 
tion takes into account of various possibilities under which the start-up singularities associated 
with parabolic problems arise. Analysis used in determining the nonuniform graded partition 
points in [5] is distinct from the one used in [10] and an example is provided in [5], which demon- 
strates that the method of Kaneko, Bey and Hou gives more sparse time partition points than 
the ones given in [10]. 
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